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Abstract-  Low  back  pain  is  a  very  common  problem  in  the 
industrialised  countries  and  its  associated  cost  is  enormous. 
Diagnosis  of  the  underlying  causes  can  be  extremely  difficult. 
Many  studies  have  focused  on  mechanical  disorders  of  the  spine. 
Digital  videofluoroscopy  (DVF)  was  widely  used  to  obtain 
images  for  motion  studies.  This  can  provide  motion  sequences  of 
the  lumbar  spine,  but  the  images  obtained  often  suffer  due  to 
noise,  exacerbated  by  the  very  low  radiation  dosage.  Thus 
determining  vertebrae  position  within  the  image  sequence 
presents  a  considerable  challenge. 

In  this  paper,  we  show  how  our  new  approach  can 
automatically  detect  the  positions  and  borders  of  vertebrae 
concurrently,  relieving  many  of  the  problems  experienced  in 
other  approaches.  First,  we  use  phase  congruency  to  relieve 
difficulty  associated  with  threshold  selection  in  edge  detection  of 
the  illumination  variant  DVF  images.  Then,  our  new  Hough 
transform  approach  is  applied  to  determine  the  moving 
vertebrae,  concurrently.  We  include  optimisation  via  a  genetic 
algorithm  as  without  it  the  extraction  of  moving  multiple 
vertebrae  is  computationally  daunting.  Our  results  show  that 
this  new  approach  can  indeed  provide  extractions  of  position 
and  rotation  which  appear  to  be  of  sufficient  quality  to  aid 
therapy  and  diagnosis  of  spinal  disorders. 

Keywords  -  Low  back  pain,  DVF,  phase  congruency,  genetic 
algorithm,  Hough  transform 

I.  INTRODUCTION 

A.  Low  Back  Pain 

The  spine  constitutes  the  central  axis  of  the  human  body 
and  can  be  divided  into  four  parts:  the  cervical,  thoracic  and 
lumbar  spine  and  the  sacroiliac  region.  The  lumbar  spine  is 
designed  to  bear  considerable  loads  and  provides  truncal 
mobility.  It  is  the  main  area  where  low  back  pain  occurs. 
However,  the  understanding  of  the  low  back  pain  is  limited 
by  the  structural  complexity  of  the  lumbar  spine  and  the 
difficulty  of  in  vivo  experiments  because  the  lumbar  spine  is 
particularly  difficult  to  access. 

Low  back  pain  and  its  associated  disability  have  appeared 
to  escalate  with  time  despite  the  considerable  technical 
advances  in  diagnosis,  treatment  and  rehabilitation.  The  cost 
of  low  back  pain  is  enormous.  For  example,  it  has  been 
estimated  that  chronic  low  back  pain  annually  results  in 
225,000  to  300,000  lumbar  surgeries  and  an  estimated  direct 
and  indirect  medical  cost  of  $75  to  $100  billion  in  the  U.  S. 
[1].  In  the  U.  K.,  the  situation  is  similar:  it  costs  billions  of 
pounds  annually  and  more  and  more  attention  has  been  paid 
to  low  back  pain  in  the  last  25  years  mainly  because  of  this 
large  cost  [2]. 

There  are,  as  yet,  no  well-accepted  standards  to  determine 
causes  of  low  back  pain.  Many  people  consider  non-specific 
low  back  pain  to  be  caused  by  abnormal  motion  and 


consequently  many  attempts  have  been  made  to  find  ways  to 
define  the  relationships  between  motion  and  low  back  pain. 

An  important  mechanical  cause  of  low  back  pain  is  spinal 
instability  which,  indeed,  might  be  one  of  the  most  common 
causes.  It  is  estimated  that  20-30  percent  of  low  back  pain 
patients  have  spinal  instability  [3].  It  is  often  held  that 
instability  may  cause  abnormal  movement  and  thus  study  of 
spinal  movement  could  help  in  its  definition  and  may  benefit 
diagnosis  and  clinical  surgery.  In  this  study,  radiographic 
methods  have  been  widely  used  in  obtaining  data  [4,  5].  In 
addition  to  these,  some  non-radiographic  effort  has  been 
devoted  to  describing  the  relationship  between  low  back  pain 
and  motion  patterns  [6,  7]. 

Spinal  instability  has  not  been  without  controversy. 
Mechanical  disorders  can  be  described  by  joint  kinematics 
and  knowledge  of  the  forces  acting  on  the  structures 
involved.  As  the  forces  exerted  on  the  lumbar  spine  are 
difficult  to  measure  in  vivo ,  clinical  studies  of  spinal 
biomechanics  have  to  focus  primarily  on  joint  kinematics.  In 
spinal  motion  analysis,  different  parameters  are  used  to 
describe  the  kinematics  [8]. 

B.  Digital  Videofluoroscopic  Imaging 

In  spinal  data  acquisition,  routine  imaging  techniques  have 
proved  to  be  unsuitable.  Due  to  the  high  radiation  dosage, 
only  a  limited  number  of  static  images  can  be  obtained  using 
plain  X-rays,  usually  in  the  neutral  position  and  at  the 
extreme  positions  of  mobility.  Consequently,  it  is  not  possible 
to  determine  the  intermediate  states  or  to  describe  the  motion 
as  the  spine  moves  from  flexion  to  extension.  Computerised 
topography  (CT)  cannot  yield  movement  information  since  it 
requires  the  patient  to  be  as  stationary  as  possible  during 
image  acquisition  whilst  magnetic  resonance  imaging  (MRI) 
is  not  yet  sufficiently  fast  in  image  acquisition  for  motion 
analysis. 

To  overcome  these  problems,  a  DVF  imaging  technique 
was  introduced  by  Breen  and  Allen  et  al.  in  1987  [9]  and  has 
undergone  considerable  refinement.  With  this  technique,  a 
series  of  dynamic  frames  of  spinal  motion  can  be  captured 
with  a  lower  X-ray  dosage  than  that  required  for  a  single 
plain  X-ray  plate  of  the  lumbar  spine.  A  typical  DVF  image 
of  the  lumbar  spine  is  shown  in  Fig.  1.  This  was  obtained 
from  a  study  of  passive  motion  in  which  the  subject  lay  on  an 
articulated  table  and  was  moved  passively  at  a  controlled  rate. 

C.  Previous  Landmark  Locating  Methods 

In  spinal  motion  study,  it  is  essential  to  locate  the 
landmarks  which  can  be  used  to  determine  the  positions  of 
the  vertebral  bodies.  This  work  was  originally  achieved 
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manually  and  consisted  of  locating,  typically,  the  corners  of 
the  vertebrae  as  anatomical  landmarks.  However,  it  is 
difficult  to  place  markers  exactly  on  the  vertebral  comers  and 
furthermore,  repeatability  cannot  be  assured.  Panjabi  et  al 
[10]  discussed  errors  that  can  arise  when  manually  marking 
X-ray  images  of  the  spine. 


Figure  1  Lateral  DVF  image  of  the  lumbar  spine 


Several  automatic  approaches  based  on  correlation  have 
been  developed  to  overcome  these  problems.  Template 
matching  methods  [11,  12]  were  proposed  but  these  may 
suffer  when  out-of-plane  motion  is  evident  or  when  occlusion 
occurs.  Other  computer  vision  techniques  have  also  been 
applied  to  vertebrae  extraction  [13,  14],  but  require  an 
excessive  amount  of  training  data.  The  latter  method  [14]  is  a 
version  of  the  active  contour  (or  snake).  One  weakness  of  a 
snake  is  that  it  depends  on  appropriate  initialization  of 
parameters.  Moreover,  snakes  cannot  solve  the 
correspondence  problems  whereby  shapes  found  might  differ 
between  frames  because  of  occlusion  or  noise  effects. 
Consequently,  errors  may  be  generated  and  would  be 
propagated  to  the  computation  of  spine  kinematics. 

In  this  paper,  we  propose  a  method  in  which  a  genetic 
algorithm  (GA)  is  combined  with  the  Hough  transform  (HT). 
It  can  represent  non-analytical  shapes  continuously  by  using 
Fourier  descriptors  (FDs),  and  thus  can  approximate  the 
model  without  distortion.  By  using  a  GA  to  search  the  Hough 
spaces,  multiple  objects  can  be  found  (here  there  are  five 
lumbar  vertebrae)  within  the  same  frame  simultaneously  and 
false  peaks  can  be  avoided  by  considering  the  intrinsic 
relationship  between  these  vertebrae.  This  approach  was 
applied  to  several  DVF  image  frames  and  the  results  are 
encouraging  as  shown  in  below. 


II.  The  hough  transform 


The  Hough  transform  [15]  has  become  one  of  the  most 
powerful  approaches  in  computer  vision.  It  has  found 
application  in  a  wide  variety  of  problems  in  machine  vision. 
Two  comprehensive  surveys  [16,  17]  of  the  HT  give  much 
evidence  of  its  attributes. 

The  philosophy  of  the  HT  is  a  mapping  from  geometric 
features  of  edge  pixels  in  an  image  to  a  multi-dimensional 
space.  Aguado  et  al.  [18]  gave  the  most  elegant  definition  of 
the  HT  based  on  the  Principle  of  Duality.  Sklansky  [19] 
proved  that  the  HT  can  provide  a  result  equivalent  to  that 
derived  by  template  matching  but  with  less  computational 


effort.  The  HT  also  inherits  advantages  such  as  immunity  to 
noise  and  occlusion. 

The  generalised  HT  (GHT)  was  first  introduced  by  Ballard 
to  extract  arbitrary  shapes  [20].  In  the  GHT,  the  model  shape 
is  represented  by  an  R-table,  which  is  a  discrete  lookup  table 
based  on  its  edge  information.  When  the  model  is  scaled  or 
rotated,  there  can  be  problems  with  aliasing  and  rounding. 
Distortions  are  inevitable  when  working  with  discrete 
representations,  as  in  digital  computers.  This  problem  was 
overcome  in  an  adaptation  of  the  GHT  [21]  where  FDs  were 
used  to  represent  the  model  shape.  This  description  gives  a 
continuous  representation  at  multiple  scales  without  the 
aliasing  problems  of  the  R-table. 

These  methods  were  originally  designed  to  extract  single 
object,  separately.  Therefore,  for  multiple  objects,  they  could 
only  be  extracted  in  sequential  stages.  However,  doing  this 
isolates  the  intrinsic  relationships  among  these  objects  that 
can  be  very  useful  in  extraction.  For  example,  during  spine 
motion,  vertebrae  are  constrained  by  the  physical  structure.  It 
might  be  possible  to  design  a  model  to  include  all  these 
objects,  but  the  dimensions  of  the  Hough  space  could  be 
excessively  large.  As  we  know,  for  a  two-dimensional  object, 
at  least  three  parameters  (rotation  and  translation  in  x  and  y 
directions  about  the  model)  are  needed  to  describe  its  pose 
and  there  will  be  fifteen  dimensions  for  five  objects. 
Furthermore,  there  is  a  possible  range  for  each  parameter  and 
this  range  is  discretised  according  to  the  resolution  during 
forming  Hough  space.  Assuming  ten  values  for  each 
parameter  implies  a  1015  search  space.  The  actual  space  in 
practice  is  even  larger  than  this,  so  traditional  exhaustive 
searches  will  be  impractical  for  so  large  scale. 

In  an  attempt  to  conquer  this  problem,  we  have  combined 
GA  with  the  HT.  The  GA  is  a  popular  method  in  optimisation 
now  and  can  be  used  to  improve  the  search  procedure.  The 
process  of  the  GA  involves  four  major  steps:  namely, 
initialization,  reproduction,  crossover  and  mutation.  The 
principle  of  the  GA  is  the  natural  phenomenon  of  “survival  of 
the  fittest”.  Limited  by  space,  we  will  not  discuss  it  in  detail 
here  and  the  reader  is  referred  to  the  extensive  literature  on  its 
implementation  and  advantages  [22].  By  introducing  a  GA 
into  the  HT,  it  not  only  can  solve  the  problem  of  the 
complexity  of  exhaustive  searching  approaches,  but  can  also 
make  it  possible  to  improve  the  performance  by  considering 
the  relationships  among  objects  in  the  fitness  function. 

III.  Results 

Edge  information  is  a  prerequisite  to  the  implementation  of 
the  HT.  As  we  discussed  earlier,  the  quality  of  DVF  images  is 
poor.  For  edge  detection,  the  Canny  operator  is  most  often 
used  and  is  regarded  as  optimal.  However,  it  cannot  provide 
good  results  on  DVF  images  of  the  lumbar  spine,  especially 
in  the  LI  and  L5  areas.  This  might  be  caused  by  the  uneven 
brightness  and  poor  contrast  within  a  single  DVF  image. 
Furthermore,  there  is  difficulty  in  selecting  optimal 
thresholds  manually  for  a  large  database  in  order  to  obtain  the 
best  edge  results.  Here,  we  have  used  phase  congruency  [23] 
and  the  results  suggest  that  it  is  appropriate.  Different  from 
Canny  and  other  gradient-based  methods,  it  utilizes  the  fact 
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that  the  feature  points  are  perceived  at  points  in  an  image 
where  the  Fourier  components  are  maximally  in  phase.  It  is  a 
dimensionless  quantity  and  is  invariant  to  changes  in  image 
brightness  and  contrast.  Furthermore,  it  provides  an  absolute 
measure  of  the  significance  of  feature  points,  and  thus  allows 
for  a  universal  threshold  that  can  be  applied  over  wide  classes 
of  images.  As  such,  difficulty  in  selecting  a  threshold  for  a 
whole  image  with  uneven  brightness,  like  DVF  images,  can 
be  solved.  Fig.  2  shows  a  comparison  between  Canny  and 
phase  congruency  where  the  source  image  is  the  central  area 
of  Fig.  1.  Visually,  phase  congruency  can  provide  better  edge 
information  and  in  particular  there  is  large  improvement  in 
boundary  detail  in  the  region  of  LI  and  L5. 


(a)  Results  with  Canny  (b)  Phase  congruency 

Figure  2  Comparisons  between  Canny  and  Phase  Congruency 


(a)  Original  curve 


At  present  the  model  shapes  were  obtained  manually,  but 
might  be  automatically  generated  in  future.  FDs  were  then 
derived  from  the  chain  code  that  can  be  easily  obtained  from 
the  model  shape  [24].  With  more  harmonics,  the 
reconstruction  becomes  increasingly  close  to  the  model. 
However,  the  total  number  of  harmonics  should  conform  to 
the  sampling  theorem.  That  is,  the  possible  values  for  number 
of  harmonics  should  be  integers  between  one  and  the  half  of 
the  sampling  rate.  Fig  3  (a)  shows  a  closed  curve  and  the  rest 
of  Fig.  3  are  its  reconstructions  with  different  harmonics.  In 
our  approach,  sixteen  FDs  were  used  to  represent  the  model. 

The  core  of  the  HT  is  how  to  form  the  Hough  space.  In  our 
study,  only  rotation  and  translations  in  the  x  and  y  directions 
were  considered.  Thus  the  Hough  space  for  each  vertebra  is 
three-dimensional.  Each  edge  point  will  vote  in  this  array  and 
the  parameters  can  be  determined  by  locating  the  maximum 
in  this  array. 

As  discussed  earlier,  a  GA  was  used  to  extract  five 
vertebrae  concurrently  by  looking  for  the  peak  values.  For 
five  vertebrae,  there  are  fifteen  parameters  altogether  and 
each  parameter  is  represented  using  six  bits,  so  the  length  of  a 
chromosome  is  90  bits.  After  trials,  the  population  number 
was  set  to  100.  A  two-point  crossover  is  used  and  its 
probability  is  0.95.  The  mutation  probability  is  0.015.  The 
program  terminates  after  2000  generations.  In  fact,  the  fitness 
value  often  converged  within  this  value.  Initialization  is 
implemented  via  random  numbers. 

The  fitness  function  is  constructed  as  the  sum  of  the  Hough 
values  corresponding  to  five  vertebrae  and  the  penalty 
constraints  on  positions.  These  constraints  are  used  to 
maintain  the  distances  between  the  neighboring  vertebral 
centers  are  within  certain  range  according  to  the  physical 
structure  of  the  lumbar  spine.  This  can  help  to  avoid  false 
peaks  caused  by  ambiguous  edge  information  (especially  for 
extraction  of  LI  and  L5).  Fig.  4  shows  the  fitness  change 
against  generation  number  in  one  example.  After  extraction, 
we  superimpose  the  results  on  the  original  image  for  an 
intuitive  observation.  Visually,  the  results  appear  very 
promising  and  one  frame  is  shown  in  Fig.  5.  Fig.  6  shows  the 
moving  patterns  of  the  lumbar  spine  when  moving  from 
neutral  posture  to  extension  position. 

550  - T - T - T - T - T - T - T - T - T - 

500 

450  -  / 


(a)  4  FDs 


(b)  8  FDs 


(c)  16  FDs  (d)  32  FDs 

Figure  3  Curve  and  its  reconstruction  with  different  FDs. 


200  400  600  800  1000  1200  1400  1600  1800  2000 

N  o.  of  generation 

Figure  4  Fitness  change  vs.  number  of  generation. 
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Figure  5  Extraction  results  of  one  DVF  image. 


—  Frame  1 

—  Frame  2 
Frame  3 

—  Frame  4 
Frame  5 


Figure  6  Extraction  results  of  multiple  frames. 

IV.  Conclusions  and  future  work 

In  this  paper,  we  show  how  the  optimized  concurrent  HT 
can  be  used  to  locate  five  lumbar  vertebrae  simultaneously 
from  DVF  images.  The  results  are  very  promising.  A  major 
advantage  is  that  it  can  consider  the  relationships  between 
objects  and  thus  it  can  help  to  avoid  false  extraction 

In  the  current  algorithm,  only  the  relative  positions 
between  vertebrae  were  considered,  but  there  is  much 
temporal  information  in  a  motion  sequence  that  can  be  used, 
for  example,  the  motion  parameters  for  each  vertebra  will  not 
change  abruptly  because  of  the  constraints  from  ligaments 
and  muscles.  Our  future  work  will  consider  incorporating  the 
temporal  information  into  our  present  approach  to  form  a 
spatial-temporal  HT  and  the  design  of  an  appropriate  fitness 
function.  The  new  approach  will  be  particularly  attractive  for 
coping  with  medical  image  sequences  of  poor  quality. 

Based  on  the  results  from  our  study,  kinematic  parameters 
of  the  lumbar  spine  can  be  easily  obtained.  Kinematics  of  the 
lumbar  spine  will  be  another  focus  of  the  future  work  and  this 
should  provide  clinicians  with  valuable  information  for 
diagnosis  of  spinal  disorders. 

Acknowledgment 

The  authors  wish  to  thank  the  Anglo-European  College  of 
Chiropractic  (AECC),  and  particularly  M.  Kondracki,  for 
kindly  providing  the  DVF  image  sequences.  Y.  Zheng 
acknowledges  the  tenure  of  a  research  studentship  supported 
by  Faculty  of  Engineering  and  Applied  Science  at  the 
University  of  Southampton. 


References 

1.  J.  A.  Porterfield  and  C.  Derosa,  Mechanical  Low  Back  Pain: 
Perspective  in  Functional  Anatomy  (2nd  edition).  W.  B.  Saunders 
Company,  1998. 

2.  Department  of  Health,  “The  prevalence  of  back  pain  in  Great  Britain  in 
1998,”  DOH  Bulletin,  1999/18,  30  June,  1999. 

3.  M.  H.  Pope  and  M.  Panjabi,  “Biomechanical  definitions  of  spinal 
instability,”  Spine,  vol.  10,  pp.  255  -  56,  1985. 

4.  N.  G.  Ogston,  G.  J.  King,  S.  D.  Gertzbein,  M.  Tile,  A.  Kapasouri  and  J. 
D.  Rubenstein,  “Centrode  patterns  in  the  lumbar  spine  baseline  studies 
in  normal  subjects,”  Spine  vol.  1 1,  pp.  591  -  595,  1986. 

5.  M.  Kanayama,  K.  Abumi,  K.  Kaneda,  S.  Tadano  and  T.  Ukai,  “Phase 
lag  of  the  intersegmental  motion  in  flexion-extension  of  the  lumbar  and 
lumbarsacral  spine,”  Spine,  vol.  21,  pp.  1416  -  1422,  1996. 

6.  T.  Steffen,  R.  K.  Rubin,  H.  G.  Baramki,  J.  Antoniou,  D.  Marchesi  and 
M.  Aebi,  “A  new  technique  for  measuring  lumbar  segmental  motion  in 
vivo :  method,  accuracy  and  preliminary  results,”  Spine,  vol.  22  (2),  pp. 
156  -  166,1997. 

7.  M.  L.  Gatton  and  M.  J.  Pearcy,  “Kinematics  and  movement  sequencing 
during  flexion  of  the  lumbar  spine,”  Clinical  Biomechanics,  vol.  14  (6), 
pp.  376  -  83,  1999. 

8.  J.  M.  Muggleton  and  R.  Allen,  “Insight  into  the  measurement  of 
vertebral  translation  in  the  sagittal  plane,”  Medical  Engineering  & 
Physics,  vol.  20,  pp.  21  -  32,  1998. 

9.  A.  C.  Breen,  R.  Allen  and  A.  Morris,  “A  computer/X-ray  method  for 
measuring  spinal  segmental  movement:  a  feasibility  study,” 
Transactions  of  the  Pacific  Consortium  for  Chiropractic  Research, 
1987,  E3.1  -E3.7. 

10.  C.  Simonis,  R.  Allen,  and  R.  Cloke,  “Spatial  analysis  of  the  movement 
of  the  spine:  Application  of  parallel  computing  in  the  field  of  spine 
biomechanics,”  Transputer  Applications,  vol.  1,  pp.  35  -  38,  1993. 

11.  J.  M.  Muggleton  and  R.  Allen,  “Automatic  location  of  vertebrae  in 
digitized  video  fluoroscopic  images  of  the  lumbar  spine,”  Medical 
Engineering  and  Physics,  vol.  19,  pp.  77  -  89,  1997. 

12.  M.  Panjabi  and  D.  Chang  and  J.  Dvorak,  “An  analysis  of  errors  in 
kinematic  parameters  associated  with  in  vivo  functional  radiographs,” 
Spine,  vol.  17(2),  pp.  200  -  205,  1992. 

13.  P.  P.  Smyth,  C.  J.  Taylor  and  J.  E.  Adams,  “Automatic  measurement  of 
vertebral  shape  using  active  shape  models,”  Image  and  Vision 
Computing,  vol.  15,  pp.  575  -  581,  1997. 

14.  M.  Brejl  and  M.  Sonka,  “Automated  initialization  automated  design  of 
border  detection  criteria  in  edge-based  segmentation,”  The  4th  IEEE 
Southwest  Symposium  on  Image  Analysis  and  Interpretation,  Austin, 
USA,  2-4  April  2000,  pp26  -  30. 

15.  P.  V.  C.  Hough,  “A  method  and  means  for  recognizing  complex 
patterns,”  U.S.  Patent  3,069,654,  1962. 

16.  J.  Illingworth  and  J.  Kittler,  “A  survey  of  the  Hough  transform,” 
Computer  Vision,  Graphics,  and  Image  Understanding,  vol.  44,  pp.  87  - 
116,1988. 

17.  V.  F.  Leavers,  “Survey:  which  Hough  transform?”  CVGIP:  Image 
Understanding,  vol.  58(2),  pp.  250  -  264,  1993. 

18.  A.  S.  Aguado,  E.  Montiel  and  M.  S.  Nixon,  “On  the  intimate 
relationship  between  the  principle  of  duality  and  the  Hough  transform,” 
Proceedings  of  the  Royal  Society-A,  vol.  456,  pp.  503  -  526,  1998. 

19.  J.  Sklansky,  “On  the  Hough  technique  for  curve  detection,”  IEEE 
Transactions  on  Computers,  vol.  C27,  pp.  923  -  926,  1978. 

20.  D.  H.  Ballard,  “Generalizing  the  Hough  transform  to  detect  arbitrary 
shapes,”  Pattern  Recognition,  vol.  13(2),  pp.  Ill  -  122,  1981. 

21.  A.  S.  Aguado,  M.  S.  Nixon  and  M.  E.  Montiel,  “Parameterizing 
arbitrary  shapes  via  Fourier  descriptors  for  evidence-gathering 
extraction,”  Computer  Vision  and  Image  Understanding,  vol.  69,  pp. 
202  -221,1998. 

22.  D.  E.  Goldberg,  Genetic  Algorithms  in  Search,  Optimization  and 
Machine  Learning.  Addison- Wesley,  October,  1998. 

23.  P.  Kovesi,  “Image  features  from  phase  congruency,”  Vidre:  Journal  of 
Computer  Vision  Research,  vol.  1(3),  pp.  3  -  26,  1999. 

24.  F.  P.  Kuhl  and  C.  R.  Giardina,  “Elliptic  Fourier  features  of  a  closed 
contour,”  Computer  Graphics  and  Image  Processing,  vol.  18,  pp.  236  - 
258,  1982. 


